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1 Introduction. 

In a recent paper |14j . Calvi and Levenberg have developed the theory of "ad- 
missible meshes" for uniform polynomial approximation over multidimensional 
compact sets. Their theory is for C , but here we restrict our attention to ffi , 
the relevant definitions and results being easily adaptable to the general case. 
We begin with introducing some notation. 

Suppose that K C M'* is compact. We will require that K not be too small, 
that is, that it is polynomial determining^ i.e., if a polynomial p is zero on K, 
then it is identically zero. This will certainly be the case if K contains some 
open ball, as will be the case for all our examples. 

We will let denote the space of polynomials of degree at most n in d real 
variables. 

Then, a Weakly Admissible Mesh (WAM) of X is a sequence of discrete 
subsets An C K such that 

IIpIIa-<C(^„)IHU„ , VpePf , (1) 

where card(A„) and the constants C{A„) have (at most) polynomial growth in 
n, i.e., 

N := dim(P^) ^(^^'^^ < card(AO = a > (2) 

and 

C{An)^0{n^), 13 >0. (3) 

Throughout the paper, \\p\\x — ^s^x^ex \p{x)\. When {C{An)} is bounded, 
C{An) < C, we speak of an Admissible Mesh (AM). It is easy to see that 
(W)AMs satisfy the following properties (cf. [H]): 

• for afhne mappings of K the image of a WAM is also a WAM, with the 
same constant C{An) 

• any sequence of sets containing An, with polynomially growing cardinali- 
ties, is a WAM with the same constants C{An) 

• any sequence of unisolvent interpolation sets whose Lebesgue constant 
grows at most polynomially with n is a WAM, C{An) being the Lebesgue 
constant itself 

• a finite union of WAMs is a WAM for the corresponding union of compacts, 
C{An) being the maximum of the corresponding constants 

• a finite cartesian product of WAMs is a WAM for the corresponding prod- 
uct of compacts, C{An) being the product of the corresponding constants. 

As shown in 14J, such meshes are very useful for polynomial approximation 
in the max-norm on K. In fact, given the classical least-square polynomial 
approximant on An to a function / G C{K), say Cnf G P,j, we have that 

\\f~Cnf\\K< (l + C(AO(l + Vcard(A„))) mm{\\f - p\\k, p eK) , (4) 
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see |141 Thm. 1]. Moreover, Fekete points (points that maximize the Vander- 
monde determinant) extracted from a WAM have a Lebesgue constant with the 
bound 

A„<iVC(A„), (5) 

that is C{An) times the classical bound for the Lebesgue constant of true (con- 
tinuum) Fekete points; see [T31 §4.4]. Recently, a new algorithm has been pro- 
posed for the computation of Approximate Fekete Points, using only standard 
tools of numerical linear algebra such as the QR factorization of Vandermonde 
matrices, cf. [711511]. 

These results show that in computational applications it is important to 
construct WAMs with low cardinalities and slowly increasing constants C(A„). 
We recall that is always possible to construct easily an AM on compact sets that 
admit a Markov polynomial inequality, as is shown in the key result [HI Thm. 
5]. There are wide classes of such compacts, for example convex bodies, and 
more generally sets satisfying an interior cone condition, but the best known 
bounds for the cardinality of the resulting meshes grow like 0{rf'^)^ where r 
is the exponent of the Markov inequality (typically r = 2 in the real case). 
This means that, for example, for a 3-dimensional cube or ball one should work 
with O^nP) points, a number that becomes practically intractable already for 
relatively small values of n (recall that the number of points determines the 
number of rows of the Least Squares (non-sparse) matrices). But already in 
dimension two, working with 0(n*) points makes the construction of polynomial 
approximants at moderate values of n computationally rather expensive. 

The properties of WAMs listed above (concerning finite unions and prod- 
ucts) suggest an alternative: we can obtain good meshes, even on complicated 
geometries, if we are able to compute WAMs of low cardinality on standard 
compact "pieces". For example, it is immediate to get a 0{n'^) WAM with 
C{An) — ©(log"* (n)) for the d-dimensional cube, as the tensor-product of 1- 
dimensional Fekete (or Chebyshev) points. 

In this paper we introduce the notion of geometric WAM, that is a WAM 
obtained by a geometric transformation of a suitable low-cardinality discretiza- 
tion mesh on some reference compact set, like the d-dimensional cube. Such 
geometric WAMs, as well as those obtained by finite unions and products, can 
be used directly for discrete Least Squares approximation, as well as for the 
extraction of good interpolation points by means of the Approximate Fekete 
Points algorithm described in [71 130j. on compact sets with various geometries. 

In Section 2 we illustrate the idea of geometric WAMs by several examples 
in M , the disk, triangles, trapezoids, and polygons. In Section 3 we prove a 
general result on the asymptotics of approximate Fekete points extracted from 
WAMs by the greedy algorithm in [71 [30] • Finally, in Section 4 we present 
some numerical results concerning discrete Least Squares approximation and 
interpolation at Approximate Fekete Points on geometric WAMs. 

2 Geometric WAMs. 

Let K and Q be compact subsets of 'Mf, 

t . Q K a surjective map (6) 
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and {Sn} a sequence of discrete subsets of Q. A geometric WAM of if is a 
sequence 

An = t{Sn) , (7) 

where the defining properties of a WAM, (HI)-© stem from the "geometric 
structure" of K, Q, t, and Sn- To make more precise this still somewhat vague 
notion, we give some illustrative examples in M , where the reference compact 
Q is a rectangle. 

2.1 The disk. 

A geometric WAM of the disk can be immediately obtained by working with 
polar coordinates, i.e., by considering the map 

t : Q = [0, 1] X [0, 2tt]^K = {x: \\x\\2 < 1} 

t(r, 0) = (r COS0, r sin0) . (8) 

We now state and prove the following 

Proposition 1 The sequence of polar grids 

r 1 1 77r If 27rfc 1 

Sn - {{rj,<f>k)}j,k = \^ + ^cos^—,0<j <n[ X I , < fc < 2n| 

gives a WAM An = t{Sn) of the unit disk, such that C(A„) — ©(log^ (n)) and 
card(A„) = + n + 1. 

Proof. Observe that given a polynomial p e P^j, when restricted to the disk in 
polar coordinates, q{r,(j3) — p{t{r,(j))), becomes a polynomial of degree n in r 
for any fixed 0, and a trigonometric polynomial of degree n'm (j) for any fixed r. 
Recall that n + 1 Chebyshev-Lobatto points are near-optimal for 1-dimensional 
polynomial interpolation, and 2n + 1 equally spaced points are near-optimal 
for trigonometric interpolation, both having a Lebesgue constant ©(log (n)); cf. 
[TSl . Now, for every p G P„ we can write 

\p(xi,X2)\ — \q{r,(j))\ — |p(rcos0, r sin(/)| < Ci log (n) max |g(rj, (/))| 

3 

where ci is independent of (j). since the {^j} are the n + 1 Chebyshev-Lobatto 
points in [0, 1]. Further 

\q{rj , 0) I < C2 log {n) max \q{rj 

k 

where C2 is independent of since the {0fc} are 2n-|- 1 equally spaced points in 
[0,27r]. Thus 

\p{xi,X2)\ < CiC2log^ (n) max|q(rj,0fc)| , V(a;i,a;2) £ K , 

i.e., An = t{Sn) = {{rj coscpk^rj sincfik)} is a WAM for the disk with C{An) = 
0(log^ (n)). We conclude by observing that the number of distinct points in 
t{Sn) is, due to the fact that r„ ~ 0, card(S'„) — 2n ~ {n + l)(2n + 1) — 2n = 
2n^+n + l. □ 
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In Figure [H we display the polar grid Sn and the WAM An in the unit disk 
for n — 8 (153 points and 137 points, respectively). In view of the structure 
of Sn and t, the points of the geometric WAM cluster at the boundary and at 
the center of the disk. Note that this technique can also be used to construct 
geometric WAMs for annuli and even a ball in higher dimensions. However we 
will not pursue this here. 
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Figure 1: The polar grid and the corresponding geometric WAM of degree n = 8 
for the unit disk. 

2.2 Polynomial maps. 

In the case that t is a polynomial map, we can give the following general 

Proposition 2 Let {i?„} be a WAM of Q, and t : Q ^ K a surjective poly- 
nomial map of degree k, i.e., t{y) ~ (ti{y), ...,td{y)) with tj G Pj!. Then 
An = t{Sn), with Sn = -Bfen, is a WAM of K such that C(A„) = C(Bfe„). 

Proof. Observing that for every x ^ K there exists some y € Q such that 
X ~ t{y), and that for every p e the composition p{t{-)) is a polynomial of 

degree < kn, we have 

\p{x)\ = \p{t{y))\ < C{Bkn)Mt{-))\U„ = C{Bkn)\\p\\tiB,„) ■ □ 

2.2.1 Triangles. 

We begin by observing that for the square Q = [—1,1]^ several WAMs Bn 
formed by good interpolation points are known, namely 

• tensor-products of 1-dimcnsional near-optimal interpolation points, such 
as Chebyshev-Lobatto points, Gauss-Lobatto points, and others; 

• the Padua points, recently studied in [S] [TTJ [T^ [TH], that are near- 
optimal for total-degree interpolation. 
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All these WAMs have C(B„) = ©(log^ (n)), but the Padua points have 
minimal cardinality N = dim(P^j) = (n + l)(n + 2)/2, whereas the cardinality 
of tensor-product points is {n + 1)^ = dim(P„ !"„)• We recall that, given the 
1-dimensional Chebyshev-Lobatto points 

C„+i = {cos(j7r/'^),0<j<n}, (9) 

the Padua points of degree n are given by the union of two Chebyshev-Lobatto 
grids 

i?„ ^ Pad„ = iCii X C-T) U (C-T X C°f2) c C„+i X C„+2 ■ (10) 

There is a simple quadratic map from the square to any triangle with vertices 
u = (ui,U2), V = (vi,V2), w — (wi,W2), namely the Duffy transformation (cf. 

m) 

tiy) = 4 - + - 2/2) + 2 - + V2) + u, (11) 

which collapses one side of the square (here 2/2 = 1) onto a vertex of the trian- 
gle (here w). By this map, the Padua points i?2n of the square (cf. (UHl)) 
are transformed to the WAM An = t(i?2n) of the triangle, with constants 
C{A„) = 0(log^ (2n)) = O(log^ (n)). The number of distinct points in t{B2n) 
is card(A„) = card(B2„) - card(C^^^i) + 1 = dim(P2„) - (n + 1) = 2^^ + 2n. In 
Figure m we display the WAMs B2n in the square and A„ in the unit simplex 
for n = 8 (153 points and 144 points, respectively). In view of the properties 
of the Padua points, the points of the geometric WAM cluster at the sides and 
especially at the vertices of the simplex. 




Figure 2: The Padua points of degree 2n — 16 and the corresponding geometric 
WAM of degree n — 8 for the unit simplex. 



2.2.2 Polynomial trapezoids. 

We consider here bidimensional compact sets of the form 

K = {x ^ {xi,X2) ■■ a<xi<b, gi{xi) < X2 < 52(a;i)} , (12) 
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where 51,52 £ A polynomial map t of degree k = u + 1 from the square 
Q = [-1, 1]2 onto is 

ii(y) = — + — ^ > ^2(y) = ^ 2/2 + ^ • (13) 

Observe that a triangle could be treated (up to an afhne tranformation) as a 
degenerate linear trapezoid. 

By Proposition 2, the Padua points Bkn of the square are mapped to the 
WAM An = t{Bkn) of the polynomial trapezoid, with C(A„) — 0{log'^{kn)) 
and card(A„) < card(Sfc„) = diTO(P^„) = {kn + l){kn + 2)/2. 

In Figure [T2l we show the WAMs An — t(i?2ri) and A„ = t{Bin) obtained 
by mapping the Padua points onto a linear trapezoid (quadratic map) and a 
cubic trapezoid (quartic map), again for n = 8 (the numbers of points are 
231 — dim(P2o) and 561 = dim(P32), respectively). As expected, we observe 
clustering at the sides and at the vertices of the trapezoids. Notice that it 
could happen that card(j4„) < card(i?fc„), namely when the graphs of 171 and 
(72 intersect at some Xi G ti{Ckn+i) (the kn + 1 Chebyshev-Lobatto points of 
[a, 6]). 




-1 -0.8 -0.6 -0.4 -0.2 0.2 D.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



Figure 3: Examples of geometric WAMs of degree 71 = 8 for a linear and a cubic 
trapezoid. 

2.2.3 Finite unions: polygons. 

A relevant example concerning finite unions is given by polygons, which are 
very important in applications such as 2-dimensional computational geometry. 
As known, any simple (no interlaced sides) and simply connected polygon with 
m vertices can be subdivided into m — 2 triangles, and this can be done by 
fast algorithms, cf. e.g. [13 [35]. Once this rough triangulation is at hand, we 
can immediately obtain a geometric WAM A„ for the polygon by union of the 
geometric WAMs constructed by mapping the Padua points on the triangles, 
with C{An) = 0(log^ (n)) and card(A„) < (m - 2){2n'^ + 2n), in view of the 
basic properties of WAMs and the results of Section 2.2.1. The points of the 
union WAM will cluster especially at the triangles common sides and vertices. 
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A similar approach is to subdivide into linear trapezoids, where we again 
obtain by finite union a geometric WAM An for the polygon with C(A„) = 
0(log^ (n)) and card(^„) — 0{rmi?) (see Section 2.2.2). In Figure H] we show 
geometric WAMs generated by subdivision into linear trapezoids of a convex 
and a nonconvex polygon, for n = 8. The method adopted, which works for 
a wide class of polygons, is that used for the generation of algebraic cubature 
points in [^Hj, where the trapezoidal panels are obtained simply by orthogonal 
projection of the sides on a fixed reference line (observe that the points cluster 
at the sides and at the line). 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 4: Examples of geometric WAMs of degree n = 8 for a convex and a 
nonconvex polygon. 



3 Approximate Fekete Points from WAMs. 

In this section we go to the general setting of polynomial interpolation in several 
complex variables. Suppose that i^T C C is a compact polynomially determin- 
ing set. If Bn = {pi,p2, ■ ■ ■ ,Pn} is a basis for P^, where N := dim(P5^) and 
Zn C if is a discrete subset of cardinality TV, then 

vdm(Z„; Bn) dci{[p{a)]peB„,aezS) (14) 

is called the corresponding Vandermonde determinant. If the determinant is 
non-zero then we may form the so-called fundamental Lagrange interpolating 
polynomials 

_ vdm(Z„\{a}U{z};g„) 

(■a{Z) .— — — . l^lOj 

vdm(Z„;B„) 

Indeed, then for any f : K ^ C the polynomial (of degree n) 

Un{z) = J2 fia)ia{z) (16) 

interpolates / at the points of Z„, i.e., n.„(a) = /(a), a G Z„. A set of points 
Fn <Z K which maximize vdm(Z„;B„) as a function of Z„, are called Fekete 
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points of degree n for K and have the special property that ||^a||if = 1 (and 
hence Lebesgue constant bounded by N) and provide a very good (often excel- 
lent) set of interpolation points for K. However, they are typically very difficult 
to compute, even for moderate values of n. 

For any determining subset An C K (thought of as a sufficiently good 
discrete model of K) the algorithm introduced by Sommariva and Vianello in 
[50] and studied by Bos and Levenberg in ^ selects, in a simple and efficient 
manner, a subset Tn C An of Approximate Fekete Points, and hence provides 
a practical alternative to true Fekete points, Fn- The optimization problem is 
nonlinear, and large-scale already for moderate values of n, but the algorithm 
is able to give an approximate solution using only standard tools of numerical 
linear algebra. 

We sketch here the algorithm, in a Matlab-like notation. The goal is extract- 
ing a maximum volume square submatrix from the rectangular N x card(A„) 
Vandermonde matrix 

V = b(a)]pes„,„e^„ , (17) 

where the polynomial basis and the array of points have been (arbitrarily) or- 
dered. The core is given by the following iteration: 

Algorithm greedy (max volume submatrix of a matrix V £ R^^*^, M > N) 

• ind = [ ] ; 

• fork = l,...,N 

— "select the largest norm column co^ij.(l^)"; ind = [ind,ik]', 

— "remove from every column of V its orthogonal projection onto coZi^ ; 

end; 

which works when V is full rank and gives an approximate solution to the NP- 
hard maximum volume problem; cf. [15] . Then, we can extract the Approximate 
Fekete Points 

= An{ind) = (A„(ii), . . . , An{iN)) ■ 

The algorithm can be conveniently implemented by the well-known QR fac- 
torization with column pivoting, originally proposed by Businger and Golub 
in 1965 [TU], and used for example by the Matlab "mldivide" or "\" operator 
in the solution of underdetermined linear systems (via the LAPACK routine 
DGEQP3, cf. 1221 ES])- 

Some remarks on the polynomial basis Bn are in order. First note that if 
B'n {qi, • ■ ■ , Qn} is some other basis of P„ then there is a transition matrix 
Tn G C ^ so that B'n = BnTN- It is easy to verify that then 

YAm{Zn;B'n) = det(r7v)vdm(2'„; S„). 

Hence true Fekete points do not depend on the basis used and also the Lagrange 
polynomials ta of (fTS]) are independent of the basis. Moreover, if An and _B„ 
are two point sets for which 

|vdm(A„,S„)| > C|vdm(B„,fi„)| 
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for some constant C, then the same inequahty holds usmg the basis B'„, i.e., 
|vdm(A„,0| >C|vdm(B„,S;)| 

since both sides scale by the same factor. 

The greedy Algorithm described above is in general affected by the basis. 
But this not withstanding, in the theorem below we show that if the initial 
set An C iiT is a WAM, then the so selected approximate Fekete points, using 
any basis, and the true Fekete points for K both have the same asymptotic 
distribution. 

Theorem 1 Suppose that K C is compact, non-pluripolar, polynomially 
convex and regular (in the sense of Pluripotential theory) and that for n — 

1,2,--- , An C K is a WAM. Let b^"\ • • • ,^>i^^} be the Approximate 

Fekete Points selected from An by the greedy Algorithm of \3U^ . using any basis 
Bn- We denote by nin the sum of the degree of the N monomials of degree at 
most n, i.e., nin = dnN/(d+ 1). Then 

(1) lim |vdm|(6S"\ b'r^)!^/™" = t{K), the transfinite diameter of K 

and 

(2) the discrete probability measures fin '■= X^jLi converge weak-* to 
the pluripotential-theoretic equilibrium measure d^K of K . 

Remark. For K = [—1,1], dfi[_i i] = ^-^^=dx; for K the unit circle S^, 

djjLgi = -^dO. If X C M'' C is compact, then K is automatically polynomi- 
ally convex. We refer the reader to [21j for other examples and more on complex 
pluripotential theory. 

By 

|vdm(4"\...,4'))| 

we will mean the Vandermonde determinant computed using the standard mono- 
mial basis. 

Note also that a set of true Fekete points Fn is also a WAM and hence we 
may take An — Fn, va which case the algorithm will select i?„ = F„ (there is 
no other choice) and so the true Fekete points must necessarily also have these 
two properties. 

Proof. We suppose that Fn = /^"\ • • ■ , /^'^} C /s: is a set of true Fekete 

points for K. Suppose further that {t^"\t2"\ ' ' ' i*jv''} C ^„ is a set of true 
Fekete points of degree n for An and that are the corresponding Lagrange 
polynomials. Then, 

\\t\\K<C{An)\\it\\A^^C{An). 

It follows that the associated Lebesgue constants 

N 

An ■■= maxE ^ ^C'(A„) 

^ i=l 

and hence, since C(v4„) is of polynomial growth, 

hm AV" = 1. 
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By Theorem 4.1 of [3] 

lim |vdm(4"\...,*i;'^)r/'"" =r(if). (18) 

n^oo 

By Zaharjuta's famous result [36], we also have 

lim |vdm(/("\...,/(;'))|i/™" (19) 

n — *oo 

Further, by [T5], we have (cf. the remarks on bases preceeding the statement of 
the Theorem) 



|vdm(5l"',...,b(;))|>^|vdm(t("\...,4"))| 



and hence 



|vdm(/("\/("\...,/(;^))| > |vdm(6("\6("\...,6(;))| 



> ^|vdm(t("\4"\...,4"))| 



Thus by and ([H]) we have 

hm \vdmib[-\bi-\.--,b^-^)\'/^-^T{K) 

n — 'oo 

as 

hm (iV!)!/"" = 1. 

n — >oo 

The final statement, that /i„ converges weak-* to dfiK then follows by the 
main result of Bcrman and Bouksom [T] (see also [1]). □ 



4 Numerical results. 

In this section, we present a suite of numerical tests concerning discrete least 
squares approximation on geometric WAMs and polynomial interpolation at 
Approximate Fekete Points extracted from them. The tests concern the 2- 
dimensional compact sets discussed in Section 2, that are the unit disk, the unit 
simplex, a linear and a cubic trapezoid, a convex and a nonconvex polygon; 
see Figures [iJIll All the tests have been done in Matlab (see ^25,), by an Intel- 
Centrino Duo T-2400 processor with 1 Gb RAM. 

In order to compute the Approximate Fekete Points, we have actually used a 
refined version of the greedy algorithm of the previous section, which is sketched 
below. 

Algorithm greedy with iterative QR refinement of the basis 

• take the Vandermonde matrix V in ((T7]); 

• K, = y* ; To = / ; 

• for fc = 0, . . . , s — 1 

Vfc = QkRk ; Uk = inv(i?fc) ; 
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Vfc+1 — VkUk ; Tk+i — TkUk ; 
end ; 

• 6=(1,...,1)*; (the choice of b is irrelevant in practice) 

• w = V*\b ; ind — find(ti; ^ 0) ; J-n — An{ind) ; 

The greedy algorithm is implemented directly by the last row above (in Mat- 
lab), irrespectively of the actual value of the vector b, and produces a set of 
Approximate Fekete Points Tn- The for loop above implements a change of 
polynomial basis from (pi, . . . ,pjv) to the nearly-orthogonal basis {qi, . . . , q^) = 
(pi, . . . ,pn)Ts with respect to the discrete inner product (/, g) — J2aGA f(^) 9i^)^ 
whose main aim is to cope possible numerical rank-deficiency and severe ill- 
conditioning arising with nonorthogonal bases (usually s = 1 or s = 2 iterations 
suffice); for a complete discussion of this algorithm we refer the reader to [30] , 




-1 -0.5 0.5 1 0.2 0.4 0.6 0.8 



Figure 5: The 45 Approximate Fekete Points of degree n = 8 extracted from 
the geometric WAMs for the disk and the simplex. 

In Figures [S][7] we display the Approximate Fekete Points of degree n — 8 
extracted form the geometric WAMs of Figures [331 The computational advan- 
tage of working with a Weakly Admissible Mesh instead of an Admissible Mesh 
is shown in Table 1, where we show the cardinalities of the relevant discrete 
sets in the case of the disk. We recall that, following [TT, Thm.5], it is always 
possible to construct an AM for a real c?-dimensional compact set that admits a 
Markov polynomial inequality with exponent r, by intersection with a uniform 
grid of stepsize 0{n~'^'^). In particular, convex compact sets have r — 2, and 
it is easily seen from the Markov polynomial inequality ||Vp(x)||2 < ?t-^||p||oo 
valid for every p e P„ and x in the disk, and the proof of [TH Thm.5], that it 
is sufficient to take a stepsize h < l/n?. The cardinality of the corresponding 
AM is then 0{n^), to be compared with the Oln?') cardinality of the geometric 
WAM (see Section 2.1). For example, at degree n = 30 an AM for the disk 
has more than 2 millions points, whereas the geometric WAM has less than 2 
thousands points. 
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Table 1: Cardinalities of different point sets in the unit disk (AFP = Approxi- 
mate Fekete Points). 

points n — b n — \Q n — \b n — 2Q n — 2b n = 30 

AM 2032 31700 159692 503868 1229072 2547452 

WAM 60 220 480 840 1300 1860 

AFP 21 66 136 231 351 496 



In Table 2, we show the Lebesgue constants of the Approximate Fekete Points 
extracted from the geometric WAMs for the disk at a sequence of degrees, using 
the algorithm above with different polynomial bases. Such Lebesgue constants 
have been evaluated numerically on a suitable control mesh, much finer than the 
extraction mesh. Without refinement iterations, the best results are obtained 
with the Logan-Shepp basis, which, as is well known, is orthogonal for standard 
Lebesgue measure (cf. [23]). On the contrary, with the monomial basis we face 
severe ill-conditioning and even numerical rank deficiency of the Vandermonde 
matrix, and we get the worst Lebesgue constants. After two refinement itera- 
tions, however, we are working in practice with a discrete orthogonal basis, the 
corresponding Vandermonde matrices are not ill-conditioned, and the Lebesgue 
constants improve and stabilize. Observe that their growth is much slower than 
that of the theoretical bound 

Table 2: Numerically evaluated Lebesgue constants (nearest integer) of the 
Approximate Fekete Points extracted from the geometric WAM in the unit 
disk, with different bases (Mon = monomial, Che — product Chebyshev, LoS 
= Logan-Shepp; in parentheses the number of refinement iterations); * means 
that the rectangular Vandermonde matrix ()17|) is numerically rank-deficient. 



basis 


n = 5 


n = 10 


n = 15 


n = 2Q 


n = 25 


n = 30 


Mon(O) 


7 


21 


34 


869 




* 


Mon(2) 


5 


24 


32 


42 


60 


81 


Che(O) 


9 


23 


30 


91 


1321 




Che(2) 


5 


24 


32 


42 


60 


81 


LoS(O) 


7 


20 


32 


52 


87 


119 


LoS(2) 


5 


24 


32 


42 


60 


81 



In Table 3 we compare, for the three test functions below, the errors (in 
the uniform norm) of discrete least squares approximation on the AM and on 
the WAM of the disk, and of interpolation at the Approximate Fekete Points 
extracted from the WAM (with two refinement iterations). The test functions 
exhibit different regularity: the first is analytic entire, the second is analytic 
nonentire (a bivariate version of the classical Runge function), the third is 
but has a singularity of the second derivatives at the origin. 

• test function 1: f{xi,X2) =cos(a;i + X2) 

• test function 2: /(xi, 2:2) = 1/(1 + 16(a;^ -I- Xj)) 
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Figure 6: The 45 Approximate Fekete Points of degree n = 8 extracted from 
the geometric WAMs for a Hnear and a cubic trapezoid. 

• test function 3: f{xi,X2) = {xf + x"^)"^^"^ 

The Vandermonde matrices have been constructed using the product Chebyshev 
basis. Both the least squares and the interpolation polynomial coefficients have 
been computed by the standard Matlab "backslash" solver, and the errors have 
been evaluated on a suitable control mesh. 



Table 3: Uniform errors of polynomial approximations on different point sets in 
the unit disk, for the three test functions above; * means computational failure 
due to large dimension (see Table 1). 





points 


n = 5 


n = 10 


n = 15 


n = 20 


n = 25 


n = 30 


test 1 


LS AM 


9E-4 


3E-10 








* 




LS WAM 


5E-4 


lE-10 


3E-15 


7E-15 


6E-15 


2E-14 




interp AFP 


lE-3 


3E-10 


2E-15 


2E-15 


2E-15 


3E-15 


test 2 


LS AM 


4E-1 


lE-1 




* 


* 






LS WAM 


5E-1 


7E-2 


5E-2 


6E-3 


4E-3 


5E-4 




interp AFP 


5E-1 


7E-2 


5E-2 


6E-3 


4E-3 


5E-4 


test 3 


LS AM 


2E-2 


2E-3 






* 






LS WAM 


2E-2 


lE-3 


7E-4 


lE-4 


2E-4 


4E-5 




interp AFP 


2E-2 


lE-3 


7E-4 


lE-4 


2E-4 


4E-5 



Notice that with the AM we have computational failure ( "out of memory" ) 
in our computing system already at degree n = 15, due to the large cardinality of 
the discrete set; see Table 1. The least squares error on the WAM is close to that 
on the AM (when comparable), which shows that geometric WAMs are a good 
choice for polynomial approximation, with a low computational cost. It is worth 
observing that in the theoretical estimate dH) we even have C (An) \/caid{ An) = 
C(n2) for the AM, and C(A„) v'card(A„) = O{nlog^ (n)) for the WAM, but 
we recall that these are overestimates, the term -\/card(A„) being in some way 
"artificial" (cf. fT^, Thm.2]). 
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Figure 7: The 45 Approximate Fekete Points of degree n = 8 extracted from 
the geometric WAMs for a convex and a nonconvex polygon. 

The following tables are devoted to numerical tests on the other domains. 
First, in Table 4 we show the Lebesgue constants of the Approximate Fekete 
Points extracted from the geometric WAMs described in Section 2. Again, the 
growth is much slower than that of the theoretical bound 

As for the simplex, the Lebesgue constant of our Approximate Fekete Points 
points is larger than that of the best points known in the literature. The case 
of the simplex has been widely studied and several specialized approaches have 
been proposed for the computation of Fekete or other good interpolation points, 
due to the relevance in the numerical treatment of PDEs by spectral-element 
and high order finite-element methods: see e.g. [20l [27l [33l [35] and references 
therein. 

On the other hand, our method for computing Approximate Fekete Points 
via geometric WAMs is quite general and flexible, since it allows to work on a 
wide class of compact sets and, differently from other computational approaches, 
up to reasonably high interpolation degrees. The good quality of the discrete 
sets used for least squares approximation and for interpolation is evidenced by 
Tables 5-9. We observe that the singularity of the third test function is in the 
interior of the domains, apart from the simplex where it is located at a vertex 
(where the discrete points cluster). This explains the better results with the 
simplex for this function. In the case of the two polygons, a change of variables 
is made in order to put the problem in the reference square [—1,1]^. 

The availability of good interpolation points in compact sets with various 
geometries has a number of potential applications. One for example is connected 
to numerical cubature. Indeed, when the moments of the underlying polynomial 
basis are known [551 EI], cubature weights associated to the Approximate Fekete 
Points can be computed as a by-product of the algorithm, simply by using the 
moments vector as right-hand side b. This gives an algebraic cubature formula, 
that can be used directly, or as a starting point towards the computation of a 
minimal formula, by the method developed in [3¥ . 

Another relevant application concerns the numerical treatment of PDEs, 
where a renewed interest is arising in global polynomial methods, such as collo- 
cation and discrete least squares methods, over general domains (see, e.g., 24.)- 
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Recently, Approximate Fekete Points have been successfully used for discrete 
least squares discretization of elliptic equations [37]. Moreover, again in the 
context of numerical PDEs, Approximate Fekete Points for polygons could play 
a role in connection with discretization methods on polygonal (non simplicial) 
meshes (see, e.g., [32] and references therein). 



Table 4: Numerically evaluated Lebesgue constants (nearest integer) of the 
Approximate Fekete Points extracted from the geometric WAMs, on different 
compact sets (product Chebyshev basis with 2 refinement iterations). 



set 


n = 5 


n 10 


n = 15 


n = 20 


n 25 


n = 30 


disk 


5 


24 


32 


42 


60 


81 


simplex 


5 


15 


25 


48 


62 


80 


linear trap 


6 


19 


34 


37 


31 


54 


cubic trap 


6 


16 


35 


45 


40 


75 


conv polyg 


7 


13 


22 


53 


53 


66 


nonconv polyg 


5 


18 


36 


35 


45 


80 



Table 5: Uniform errors of polynomial approximations in the unit simplex, for 
the three test functions above. 





points 


71 = 5 


71 = 10 


77 = 15 


71 = 20 


77^ = 25 


77^ = 30 


test 1 


LS WAM 


7E-7 


8E-15 


3E-15 


4E-15 


4E-15 


6E-15 




interp AFP 


2E-6 


2E-14 


lE-15 


3E-15 


3E-15 


5E-15 


test 2 


LS WAM 


2E-2 


5E-4 


lE-5 


4E-7 


lE-8 


4E-10 




interp AFP 


5E-2 


2E-3 


4E-5 


2E-6 


3E-8 


2E-9 


test 3 


LS WAM 


7E-4 


5E-6 


4E-7 


8E-8 


2E-8 


7E-9 




interp AFP 


8E-4 


2E-5 


lE-6 


2E-7 


6E-8 


3E-8 



Table 6: As in Table 5, for the linear trapezoid of Figure 3. 





points 


71 = 5 


n = 10 


77 = 15 


71 = 20 


77^ = 25 


71 30 


test 1 


LS WAM 


3E-3 


5E-9 


lE-13 


3E-15 


4E-15 


9E-15 




interp AFP 


8E-3 


2E-8 


3E-13 


4E-15 


3E-15 


4E-15 


test 2 


LS WAM 


2E-1 


2E-1 


lE-1 


3E-2 


lE-2 


5E-3 




interp AFP 


3E-1 


2E-1 


2E-1 


3E-2 


2E-1 


lE-2 


test 3 


LS WAM 


3E-2 


4E-3 


2E-3 


5E-4 


2E-4 


lE-4 




interp AFP 


5E-2 


4E-3 


3E-3 


5E-4 


3E-4 


2E-4 
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Table 7: As in Table 5, for the cubic trapezoid of Figure 3. 





points 


n — 5 


71 = 10 


71 = 15 


n = 20 


71 = 25 


71 30 


Leo b J_ 


LS WAM 


2E-3 




6E-14 


3E-15 


4E-15 


5E-15 




interp AFP 


6E-3 


lE-8 


lE-13 


5E-15 


3E-15 


4E-15 


test 2 


LS WAM 


4E-1 


2E-1 


6E-2 


3E-2 


9E-3 


5E-3 




interp AFP 


5E-1 


2E-1 


7E-2 


5E-2 


lE-2 


6E-3 


test 3 


LS WAM 


3E-2 


3E-3 


9E-4 


5E-4 


2E-4 


2E-4 




interp AFP 


6E-2 


5E-3 


9E-4 


7E-4 


2E-4 


2E-4 


Table 8: 


: As in table 5 for the 


convex 


polygon of Figure 


4 and the three test 


functions f{2xi — 1,2x2 — 1) above. 












points 


n — 5 


71 = 10 


71 = 15 


71 = 20 


71 = 25 


71 = 30 


test 1 


LS WAM 


7E-4 


lE-9 


7E-15 


9E-15 


lE-14 


2E-14 




interp AFP 


lE-3 


4E-9 


6E-15 


4E-15 


4E-15 


5E-15 


test 2 


LS WAM 


4E-1 


lE-1 


4E-2 


2E-2 


4E-3 


lE-3 




interp AFP 


5E-1 


lE-1 


4E-2 


2E-2 


9E-3 


3E-3 


test 3 


LS WAM 


2E-2 


2E-3 


6E-4 


3E-4 


lE-4 


9E-5 




interp AFP 


2E-2 


2E-3 


6E-4 


3E-4 


lE-4 


8E-5 




Table 9: As in 


Table 8 


for the ] 


Qonconvex 


polygon 


of Figure 


4. 




points 


n = 5 


n = 10 


71 = 15 


71 ^ 20 


71 = 25 


71 = 30 


test 1 


LS WAM 


5E-4 


3E-10 


lE-14 


2E-14 


3E-14 


4E-13 




interp AFP 


6E-4 


5E-10 


3E-15 


3E-15 


3E-15 


4E-15 


test 2 


LS WAM 


4E-1 


2E-1 


5E-2 


2E-2 


5E-3 


lE-3 




interp AFP 


6E-1 


2E-1 


5E-2 


2E-2 


5E-3 


2E-3 


test 3 


LS WAM 


2E-2 


3E-3 


7E-4 


3E-4 


lE-4 


9E-5 




interp AFP 


4E-2 


3E-3 


8E-4 


3E-4 


lE-4 


7E-5 
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